Cite this Article
Song Chi, Kong Xiang-Shan, Liu C S. First-principles studies on carbon diffusion in tungsten*

Project supported by the National Key Research and Development Program of China (Grant No. 2018YFE0308102), the National Natural Science Foundation of China (Grant Nos. 11735015 and 51771185), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant No. 17KJB140008), and Jinling Institute of Technology, China (Grant Nos. jit-fhxm-201601 and jit-b-201616).

. Chinese Physics B, 2019, 28(11): 116106  
Permissions
First-principles studies on carbon diffusion in tungsten*

Project supported by the National Key Research and Development Program of China (Grant No. 2018YFE0308102), the National Natural Science Foundation of China (Grant Nos. 11735015 and 51771185), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant No. 17KJB140008), and Jinling Institute of Technology, China (Grant Nos. jit-fhxm-201601 and jit-b-201616).

Song Chi1, Kong Xiang-Shan2, †, Liu C S2
College of Science, Jinling Institute of Technology, Nanjing 211169, China
Key Laboratory of Materials Physics, Institute of Solid State Physics, Chinese Academy of Sciences, Hefei 230031, China

 

† Corresponding author. E-mail: xskong@issp.ac.cn

Project supported by the National Key Research and Development Program of China (Grant No. 2018YFE0308102), the National Natural Science Foundation of China (Grant Nos. 11735015 and 51771185), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant No. 17KJB140008), and Jinling Institute of Technology, China (Grant Nos. jit-fhxm-201601 and jit-b-201616).

Abstract

The carbon diffusivity in tungsten is one fundamental and essential factor in the application of tungsten as plasma-facing materials for fusion reactors and substrates for diamond growth. However, data on this are quite scarce and largely scattered. We perform a series of first-principles calculations to predict the diffusion parameters of carbon in tungsten, and evaluate the effect of temperature on them by introducing lattice expansion and phonon vibration. The carbon atom prefers to occupy octahedral interstitial site rather than tetrahedral interstitial site, and the minimum energy path for its diffusion goes through a tetrahedral site. The temperature has little effect on the pre-exponential factor but a marked effect on the activation energy, which linearly increases with the temperature. Our predicted results are well consistent with the experimental data obtained at high temperature (>1800 K) but significantly larger than the experimental results at low temperature (<1800 K).

PACS: 61.80.-x;61.82.Bg
Keyword:metals and alloys;diffusion;computer simulations
1. Introduction

Diffusion of carbon atoms in metals is of great scientific and technological interests because it is closely related to many important phenomena in materials science.[1–5] A typical example is that the diffusion of carbon controls kinetics of carbonization of parent metal and alloy elements.[3,4] In addition, diffusion of impurity carbon significantly affects its surfaces and grain boundary segregation properties, which in turn leads to changing mechanical properties of metals.[5,6] Recently, carbon behavior in tungsten has attracted a great deal of attention. Firstly, both tungsten and carbon have been selected as the plasma-facing materials (PFM) in the international thermonuclear experimental reactor (ITER).[7,8] During operation of this device, carbon will be inevitably introduced in tungsten and thus become one of the most frequent impurity atoms in it. Secondly, the drawback of tungsten as PFM in future fusion reactors is its high ductile-to-brittle transition temperature (DBTT) and therefore high brittleness at the operation temperature.[9] The introduce of fine and uniform dispersion of carbides in tungsten is an effective way to improve ductility by lowering DBTT.[9,10] In addition to the application in fusion reactors, tungsten has also been considered as the best substrates for diamond growth.[11] The diffusion behavior of carbon in tungsten substrates can affect the nucleation and surface morphology of diamond. For these reasons, it is very important to know the diffusion behavior of carbon in tungsten.

The diffusion of carbon in tungsten has been thoroughly investigated experimentally for over 50 years.[12–20] However, data on this are still scarce and, more importantly, the values of diffusion parameters obtained are different. As shown in Table 1, most of the pre-exponential factors are in the magnitude of 10−7 m2/s, showing a good consistency. While the activation energy varies from 1.64 eV to 2.56 eV for different measurement temperature ranges. At present, it is very hard to give a reliable value for carbon diffusivity in tungsten, and data for this is urgently required. Fortunately, it is now possible to calculate the diffusivity of light elements in metals by the first-principles density functional theory coupled with harmonic transition state theory.[21] The accuracy of these calculations are close to and sometimes better than that available from experiments. Many successful cases have been reported in the past decade, such as hydrogen and oxygen in nickel[22] and iron,[23,24] carbon in palladium and palladium alloys,[25] hydrogen in tungsten,[26–29] and so on. Recently, several researchers have used the first principles calculations to investigate the diffusion path of carbon in tungsten.[30–33] Three different activation energies of carbon diffusion were given, i.e., 1.24 eV,[30] 1.46 eV,[31] and 1.47 eV.[32,33] Moreover, Liu et al.[31] have further predicted the carbon diffusivity as D = 2.56 × 10−7 exp (−1.46 eV/kT) according to the theory of Zener and Wert (here k is Boltzmann’s constant and T is the temperature). Unfortunately, there is no comparison with experimental data in these studies. As shown in Table 1, the predicted pre-exponential factor by Liu et al.[31] is comparable with experimental results. However, all predicted activation energies are much smaller than the experimental data. This large difference may result from the temperature effect, which has been demonstrated to be an important factor in the accurate calculation of the diffusion parameters by first-principles computations.[22–25,28,29] Therefore, in this paper, we are going to carry out a systematic first-principles energy and vibration spectrum calculations to investigate the temperature effect on the carbon diffusion in bulk bcc tungsten.

Table 1.

The pre-exponential factors and activation energies obtained by experimental measurements and first-principles calculations

.
2. Computational models and methods
2.1. Computational models

The diffusivity of interstitial solute atoms in solid solutions can be expressed as[34]

where λ and Γ are the jump length and rate between adjacent sites of the diffusing particle, respectively. The jump rate Γ can be calculated by two distinct methods.

One method is based on Eyring’s theory of the activated complex, involving time-consuming calculations of vibrational frequencies,[35–38] which is denoted by hTST. In this method, the jump rate (denoted by ΓhTST) is given by

where h denote Planck’s constant, and ZTS and ZGS are partition functions for the transition and ground states, respectively. Their ratio can be written as
within the framework of the harmonic approximation. Here N is the number of the vibrating atoms, denotes the vibrational frequencies in the transition state (TS), and corresponds to the frequencies in the ground state (GS). Eact is the activation energy, defined as the energy difference between the transition and ground states, i.e.,
The other method is based on the diffusion theory proposed by Wert and Zener (denoted by WZ),[39] where the jump rate ΓWZ can be approximately estimated by
where n is the number of the nearest-neighboring interstitial positions, and m is the mass of the interstitial solute atom.

The temperature effect is taken into account by considering the thermal expansion and phonon vibration. Lu et al.[40] have proposed an equation to describe the lattice thermal expansion of tungsten, whose calculated results are fully consistent with the experimental data reported. Therefore, this equation can be used to correlate the temperature with the lattice parameter value. Similarly to previous works,[25,28,29] the relevant lattice constant (aT) at temperature T can be calculated by multiplying the DFT-optimized lattice constant a0 at 0 K by the ratio obtained using the equation of Lu et al., i.e.,

In the quasi-harmonic approximation, the vibration free energy of system A with a given lattice constant aT at temperature T can be expressed as
where is the vibration modes of the system A at temperature T. The total energy of the system A can be calculated by
Here is the static energy, obtained directly from the DFT calculations.

2.2. Computational methods

The present calculations are performed by the Vienna ab-initio simulation package (VASP) code based on density functional theory (DFT). The projector augmented wave (PAW) pseudo-potentials are used to describe interactions between ions and valence electrons.[41] For the exchange–correlation potential, the generalized gradient approximation as parameterized by Perdew–Wang (PW91) is used.[42,43] A supercell containing 128 lattice points (4 × 4 × 4) is used. Throughout the calculations, the atomic position is allowed to relax while the shape and volume of the supercell are fixed. The plane wave cutoff of 500 eV and k-point grid density of 5 × 5 × 5, generated by the Monkhorst–Pack method,[44] are employed, which give a convergence of the total energy within 0.001 eV per atom. The structural optimization is truncated when forces on all atoms in the supercell are less than 0.01 eV/Å. Lattice vibrations are calculated using the harmonic approximation and the frozen phonon approach. Similar to the previous work,[28] localized vibrations of the interstitial carbon atom and its neighboring tungsten atoms are assumed to be decoupled from vibrations of other metal atoms, i.e., only the interstitial carbon atom and its neighboring tungsten atoms are allowed to vibrate.

3. Results and discussion
3.1. Verification of the calculations

The calculated lattice constant of the tungsten at 0 K is 3.177 Å, which agrees well with the previous calculation results[29–31,33] and is slightly larger than the experimental value (3.165 Å[45]). Using Eq. (6), the relevant lattice constants (aT) at temperature T are calculated and presented in Fig. 1, which is used in the following calculation. Note that the lattice expansion corresponds to the isotropic tensile strain. Therefore, the corresponding lattice strain εaT are also presented in Fig. 1.

Fig. 1. The relevant lattice constant (aT) and their corresponding isotropic tensile strain.

To verify the accuracy of our calculations, the solution and diffusion properties of interstitial carbon in tungsten at 0 K are calculated. Being small enough, carbon will occupy interstitial sites in bcc metal. Two possible interstitial sites, i.e. tetrahedral (Tet) and octahedral (Oct) sites, are considered. The formation energy is calculated by Ef = EWC – EW – EC, where EWC is the total energy of the tungsten lattice with a carbon atom at interstitial site, EW is the total energy of tungsten perfect supercell, and EC is the energy of carbon isolated. Using this equation, the formation energies of carbon in tetrahedral and octahedral sites are calculated to be −4.45 eV and −5.92 eV, respectively. This means that the carbon atom prefers to occupy octahedral interstitial site rather than tetrahedral interstitial site. For seeking the saddle point in the process of the carbon atom migration, the nudged elastic band (NEB) method is utilized. The interpolation points of the system are set 5 to construct a minimum energy path. According to our calculation, the minimum energy path for the interstitial carbon diffusion is found directly through a Tet, i.e., the path Oct→Tet→Oct shown in Fig. 2. Further calculations on the vibrational frequencies observed only one imaginary frequency at the state with a interstitial carbon atom at Tet. This confirms that the Tet is a saddle point for the interstitial carbon diffusion. Thus, the migration energy is simply given by the difference in energy between the states of interstitial carbon atom at Oct and Tet, i.e., 1.47 eV. The pre-exponential factor calculated by Eq. (5) is 2.54 × 10−7 m2/s. These calculated results are in good agreement with previously reported results (see Table 1).[31–33]

Fig. 2. The schematic diagram of possible interstitial sites for carbon in bcc tungsten: tetrahedral and octahedral interstitial sites, which are denoted by the small red and green spheres, respectively. The tungsten atoms are represented by the large sphere. Note that the tungsten atoms marked by 1 and 2 are allowed to vibrate in the frequency calculation.
3.2. The temperature effects on the diffusion

The temperature effect on the activation energy can be evaluated by substituting Eq. (8) into Eq. (4), i.e.,

which can be divided into two parts, i.e., the contribution of the static energy
and the contribution of the vibration free energy
Figure 3 presents the calculated Eact(aT), , and . It can be clearly seen that decreases slowly and linearly with the temperature and its value decreases about 0.12 eV when the temperature increases from 300 K to 2700 K. On the contrary, Eact(aT) shows a positive dependence on the temperature, somewhat analogous to the behavior of . This suggests that the major contribution to the increase of the diffusion activation energies across the temperature comes from the vibration free energy.

Fig. 3. The temperature-dependent activation energies and their corresponding static energy and vibration free-energy contribution. The green stars represent the static energy contribution calculated by Eq. (14).
Table 2.

The carbon interstitial atom induced lattice stress (in GPa) in bcc tungsten. The compressive stress is negative and the tensile stress is positive.

.

Similar to the case of hydrogen in tungsten,[28,46] the dependence of on temperature results from the lattice expansion and can also be understand by the linear elasticity theory. According to linear elasticity theory, the solution energy of interstitial solute atom in system A under isotropic tensile strain ε can be expressed as

Here V is the volume of the supercell at equilibrium, and is the average stress induced by the interstitial solute. ESolA is the solution energy of the solute in the system A without strain, which can be calculated by
where and represent the static energy of the system A with and without solute, and ESoliso is the energy of the isolated solute atom. Thus, for interstitial carbon diffusion in tungsten, can be calculated by
according to Eqs. (10), (12), and (13). Here and are the lattice stress induced by carbon located in Oct and Tet, respectively, whose values are summarized in Table 2. The calculated values of using Eq. (14) are presented in Fig. 3, which are in close agreement with the DFT data. Therefore, the static energy contribution of the activation energy shows a linear-like negative dependence on the temperature.

Lattice expansion in general leads to a softening of vibrational modes, which in turn affects the vibration free energy and the partition functions (see section 2.1). In previous works on carbon diffusion in palladium[25] and hydrogen diffusion in tungsten,[28] researchers found that the softening of vibrational modes induced by lattice expansion has little effect on the temperature-dependent behavior of the vibration free energy. Similar results are also found in our present work. As shown in Fig. 4, the values of and ZTS/ZGS calculated using and , respectively, are almost identical. Thus, in the calculation of the interstitial solute diffusivity, one can use the vibration frequencies of the system with the a0 at all temperatures in order to save substantial amounts of computational time without significant loss of accuracy. Note that the explicitly values of and ZTS/ZGS are still used in this work.

Fig. 4. The effect of the softening of vibrational modes on and ZTS/ZGS.

Using the two distinct methods of the calculation of the jump rate described in Section 2, the diffusion coefficients of carbon in tungsten with and without the temperature effect are calculated, i.e., the jump rate from the Wert–Zener model with and without the temperature effect (denoted by WZ-TE and WZ-UNTE, respectively) and the jump rate from the harmonic transition state theory of Eqs. (2) and (3) with and without the temperature effect (named by hTST-TE and hTST-UNTE, respectively). For the situation with the temperature effect, the above temperature-dependent lattice constant aT, activation energy Eact(aT) and vibration modes are used. When the temperature effect is excluded, we only use the diffusion activation energy (1.47 eV) and the vibration modes of the system with the DFT-optimized lattice parameter a0 = 3.177 Å. The calculated results of the carbon diffusion coefficients of these four cases are displayed in Fig. 5(a). It can be seen that the predicted diffusivity using these four different methods are quite close to each other. In order to more clearly evaluate their difference, the ratios of them to the case of hTST-TE are calculated and presented in Fig. 5(b). It can be seen that the calculated diffusivities with the WZ model are around twice as much as that with the hTST model. Meanwhile, both the WZ and hTST models would overestimate the diffusion coefficient at high temperatures when the temperature effect is ignored during the calculation. On the whole, the difference between them are less than an order of magnitude. This suggests that the temperature effect on the carbon diffusivity in tungsten is not significant. That is because the temperature has little effect on the pre-exponential factor (see Fig. 6). Meanwhile, although the increase in temperature leads to a change in activation energy, this change contributes less to the diffusion coefficient at high temperatures. Therefore, the predicted carbon diffusivity can be approximately fitted by the Arrhenius equation D = D0 exp(−Eact/kT) with D0 = 2.58 × 10−7 m2/s and Eact = 1.57 eV.

Fig. 5. The predicted diffusivity of carbon in tungsten with and without the temperature effect and the corresponding pre-exponential factors and activation energies.
Fig. 6. The comparison between the experimental measurements[12–19] and our first-principles calculated diffusivity (a) and (b), activation energies (c) and pre-exponential factors (d) of carbon in tungsten. The gray area indicates the range of T > 1800 K.
3.3. Comparison with experiments

The agreement with experimental results is an extremely important criterion to verify theoretical predictions. For the interstitial diffusion of light elements in metals, the diffusivity data obtained at high temperatures are usually believed to be more reliable since they are likely less influenced by both surface and trapping effects. However, it is difficult to experimentally determine the temperature range without trapping effects due to the lack of the benchmarks. A comparison of experimental and theoretical data may hopefully solve this problem. Figures 6(a) and 6(b) exhibit the available experimental data of the carbon diffusivity in tungsten collated from the literature together with our calculated results here. Particularly, the pre-exponential factors and activation energies obtained by experimental measurement (Table 1) and first-principles calculation are also plotted in Figs. 6(c) and 6(d) to make a comparison and clearly display their difference. It can be seen from Figs. 6(a) and 6(b) that our predicted diffusion coefficients are in very good agreement with the experimental data when the temperature is above 1800 K. Moreover, our predicted activation energies are in the range of 1.58–1.67 eV at temperatures 1800–2700 K and increase with the temperature, which are comparable to those obtained by fitting the experimental data of carbon diffusion coefficients at high temperature (1.64 eV and 1.76 eV in the temperature ranges 1773–2073 K and 2073–3073 K, respectively, see Fig. 6(c). Unlike activation energy, our predicted diffusion pre-factors show no significant change with the temperature, which are around 1.3 × 10−7 and 2.5 × 10−7 for the models of hTST and WZ, respectively. They are roughly equivalent to most of the experimental results (see Fig. 6(d)). These consistencies not only indicate the reliability of our theoretical predictions but also shed some light on the experimental data at temperatures above 1800 K, less affected by trapping effects. When the temperature below ∼1800 K, our predicted diffusivity is larger than the experimental values. The difference between the experimental measurements and the theoretical predictions becomes larger as the temperature decreases. Moreover, as shown in Fig. 6(c), experiments give activation energies of 2.09 eV at 1523–1723 K, 2.32 eV at 1473–1873 K, and 2.56 eV at 1073–1143 K, which are obviously larger than our theoretical predictions. This indicates that the diffusion of carbon in the temperature range of T < 1800 K could be significantly affected by trapping effects. The defect-trapping effects on carbon diffusivity are very complex and related works are in progress.

4. Conclusion

We have performed a series of first-principles calculations to predict the carbon diffusion parameters in tungsten, namely, pre-exponential factor, activation energy, and diffusivity. The temperature effects on these diffusion parameters are evaluated by considering the lattice expansion and phonon vibration. It is found that the carbon atom prefers to occupy octahedral interstitial site rather than tetrahedral interstitial site, and the minimum energy path for its diffusion goes through a tetrahedral site. The predicted pre-exponential factor is around 3 × 10−7 m2/s over the investigated temperature range of 300–2700 K, which is not affected by temperature and consistent with the experimental results. The predicted activation energies linearly increase from 1.47 eV to 1.67 eV at temperatures 300–2700 K, which are consistent with the experimental data obtained at high temperatures (>1800 K) but much lower than those obtained at low temperatures (<1800 K). The temperature effect on the carbon diffusivity in tungsten is not significant because the change of activation energies with temperature contributes less to the diffusion coefficient at high temperatures. Similarly to the activation energies, our predicted diffusivity data are very good consistent with the high-temperature experimental values (>1800 K) but larger than the low-temperature experimental results.

Reference
[1] Teschner D Borsodi J Wootsch A Revay Z Havecker M Knop-Gericke A Jackson S D Schlogl R 2008 Science 320 86
[2] Hiraoka Y Imamura K Kadokura T Yamamoto Y 2010 J. Alloys Compd. 489 42
[3] Zhang S Wang S Li W Zhu Y Chen Z 2011 J. Alloys Compd. 509 8327
[4] Song G M Zhou Y Wang Y J 2002 J. Mater. Sci. 37 3541
[5] Hiraoka Y 1990 Trans. JIM 31 861
[6] Zhou H B Jin S 2013 Chin. Phys. 22 076104
[7] Federici G Andrew P Barabaschi P et al. 2003 J. Nucl. Mater. 313-316 11
[8] Fu B Q Lai W S Yuan Y 2013 Chin. Phys. 22 126201
[9] Rieth M Dudarev S L Gonzalez de Vicente S M et al. 2013 J. Nucl. Mater. 432 482
[10] Xu M Y Luo L M Lin J S Xu Y Zan X Xu Q Tokunaga K Zhu X Y Wu Y C 2018 J. Alloys Compd. 766 784
[11] Yang T Wei Q Qi Y Yu Z 2015 Diam. Relat. Mater. 52 49
[12] Kovenskii I 1964 Diffus. Body Cent. Cubic Met. Pap. Int. Conf. Gatlinburgh TN 283
[13] Shepela A 1972 J. Less-Common. Met. 26 33
[14] Aleksandrov L N Shchelkonogov V 1964 Sov. Powder Metall. Met. Ceramic. 288
[15] Nakonechnikov A I Pavlinov L V Bykov V N 1966 Phys. Met. Metallogr. 22 73
[16] Rawlings K Foulias S Hopkins B 1981 Surf. Sci. 109 513
[17] Aleksandrov L N 1963 Int. Chem. Eng. 3 108
[18] Shchelkonogov V Y 1978 Elektron. Svoistva Tverd. Tel. Fazovye Prevrashch. 115
[19] Eckstein W Shulga V I 1999 Nucl. Instrum. Meth. 153 415
[20] Schmid K Roth J 2002 J. Nucl. Mater. 302 96
[21] Qiao L Wang S M Zhang X M Hu X Y Zeng Y Zheng W T 2014 Chin. Phys. 23 086802
[22] Nama H O Hwang I S Lee K H Kim J H 2013 Corros. Sci. 75 248
[23] Chen Q Yao Q Liu Y L Han Q F Ding F 2017 Int. J. Hydrog. Energy 42 11560
[24] Shang S L Fang H Z Wanga J Guo C P Wang Y Jablonski P D Du Y Liu Z K 2014 Corros. Sci. 83 94
[25] Ling C Sholl D S 2009 Phys. Rev. 80 214202
[26] Heinola K Ahlgren T 2010 J. Appl. Phys. 107 113531
[27] Qiu M J Zhai L Cui J C Fu B Q Li M Hou Q 2018 Chin. Phys. 27 073103
[28] Kong X S Wang S Wu X B You Y W Liu C S Fang Q F Chen J L Luo G N 2015 Acta Mater. 84 426
[29] Liu Y L Ding F Luo G N Chen C A 2017 Nucl. Fusion 57 126024
[30] Nguyen-Manh D 2009 Adv. Mater. Res. 59 253
[31] Liu Y L Zhou H B Jin S Zhang Y Lu G H 2010 J. Phys.: Condens. Matter 22 445504
[32] Kong X S You Y W Song C Fang Q F Chen J L Luo G N Liu C S 2012 J. Nucl. Mater. 430 270
[33] Becquart C S Domain C 2012 Curr. Opin. Solid. State. Mater. Sci. 16 115
[34] Philibert J 1991 Atom Movements: Diffusion and Mass Transports in Solids, Monographies de Physique Les Ulis, France éditions de Physique
[35] Eyring H 1935 J. Chem. Phys. 3 107
[36] Eyring H 1938 Trans. Faraday Soc. 34 41
[37] Wigner E 1932 Z. Phys. Chem. Abt. 19 203
[38] Wigner E 1938 Trans. Faraday Soc. 34 29
[39] Wert C Zener C 1949 Phys. Rev. 76 1169
[40] Lu X G Selleby M Sundman B 2005 CALPHAD: Comput. Coupling Phase Diagrams Thermochem. 29 68
[41] Blöchl P E 1994 Phys. Rev. 50 17953
[42] Perdew J P Chevary J A Vosko S H Jackson K A Pederson M R Singh D J Fiolhais C 1992 Phys. Rev. 46 6671
[43] Perdew J P Chevary J A Vosko S H Jackson K A Pederson M R Singh D J Fiolhais C 1993 Phys. Rev. 48 4978
[44] Monkhorst H J Pack J D 1976 Phys. Rev. 13 5188
[45] Kittel C 1996 Introduction Solid State Physics 7 New York Wiley
[46] Zhou H B Jin S Zhang Y Lu G H Liu F 2012 Phys. Rev. Lett. 109 135502